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Abstract 

A second-order numerical implementation is given for recently derived non¬ 
linear wave equations for general relativity. The Gowdy T^ cosmology is used 
as a test bed for studying the accuracy and convergence of simulations of 
one-dimensional nonlinear waves. The complete freedom in space-time slicing 
in the present formalism is exploited to compute in the Gowdy line-element. 
Second-order convergence is found by direct comparison of the results with 
either analytical solutions for polarized waves, or solutions obtained from 
Gowdy’s reduced wave equations for the more general unpolarized waves. 
Some directions for extensions are discussed. 
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I. INTRODUCTION 


The gravitational wave observatories LIGO and VIRGO |TT|Jl[| presently under construc¬ 
tion have given added impetus towards accurate prediction of gravitational wave forms from 
binary coalescence of neutron stars and black holes. The late stages of the spiral infall are 
predominantly targeted through large scale simulations in numerical relativity. The gravita¬ 
tional wave interactions have motivated a description of gravity by nonlinear wave equations 
in a tetrad approach ( ^j). In this formulation, the strictly hyperbolic nature of the equa¬ 
tions is independent of the particular choice of foliation. The foliation is governed by four 
tetrad lapse functions which are algebraically related to the more familiar Hamiltonian lapse 
and shift functions. 

In this paper, a numerical implementation of the nonlinear wave equations is given, and 
the performance is studied in the Gowdy T^ cosmology. The Gowdy T^ is a good test 
bed for one-dimensional nonlinear wave motion of both polarized and unpolarized waves. 
The numerical scheme is one-dimensional, second-order in time and spectrally accurate in 
space. The SO(3,1, R)-connections are evolved by implementation of the four-divergence 
of the Riemann tensor and the Lorentz gauge on the connections. The tetrad elements are 
evolved by the equations of structure, and treated as a system of ordinary differential equa¬ 
tions in which the fundamental matrix is a finite Lorentz transformation. The Gowdy T^ 
line-element has been incorporated in the numerical scheme, so as to enable direct compar¬ 
ison of the computed solution with either analytical or numerical solutions obtained from 
Gowdy’s reduced wave equations. Gonvergence results are presented for both polarized and 
unpolarized Gowdy waves. 

The nonlinear wave equations for the connections, of the tetrad elements, {(e^)a} 
satisfy nonlinear wave equations of the Yang-Mills type. They follow from a Lorentz gauge 
on the connections. In vacuo, they simplify to 

- [o;'", = 0. (1) 

Here, □ = V'^Vc with Va the SO(3,l,R) -gauge covariant derivative satisfying = 


2 


Va(e^)fe + = 0. The tetrad elements define the metric by gab = Here, 

contraction over the Greek indices is through = diag(—1,1,1,1). Initial data in this 
second-order formalism includes the initial data for the familiar first-order Hamiltonian 
formalism. Additional initial data are for example initial values for the tetrad; with sufficient 
smoothness these data generate initial data for the connections and the Riemann tensor. A 
procedure for obtaining smooth initial tetrad data from the metric is included. 

In Section 2, the Gowdy waves and their reduced wave equations are discussed. Section 
3 presents the integration scheme, together with the procedure for computing initial data 
for the tetrad elements. Section 4 discusses the gauge conditions in the Gowdy line-element. 
Discussion of the simulations and conclusions are given in Section 5. 

II. GOWDY WAVES 

Gowdy cosmologies describe an extensively studied class of universes with compact space¬ 
like hypersurfaces with two Killing vectors, and ds- By choice of boundary conditions, 
the space-like hypersurfaces are homeomorphic to either the three-torus, the three-handle 
or the three-sphere. The three-torus describes a semi-infinite time-evolution of universes 
collapsing into, or beginning with a singularity. 

The Gowdy three-torus cosmology has recently been investigated numerically by Berger 
and Moncrief P|, and can be given in the line-element 

ds2 = e-tei(-e-2-dr2 + d6»2) + dE2, (2) 

where 

dS^ = e“’'[e'^dcr^ -|- 2e^Qdad6 + + e~^)d5‘^]. (3) 

Here A = X{t,6),P = P{t,6) and Q = Q{t,9), by invariance in the angular coordinates a 
and 6. 

In small amplitude waves, e~'^P and e~'^Q are the amplitudes of the e+ and gravita¬ 
tional wave polarization tensors (Eqn. 3.17 in ^j). The quantities P and Q satisfy Gowdy’s 
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reduced wave equations, which may be given as 0: 


Q,, =e-^^Qee-2{PrQr-e-^^PgQg), 

P,, =e-^^Pgg + e-^^{Ql-e-^^Ql). 

The 0—momentum and Hamiltonian constraints, respectively, are given by 


( 4 ) 


Ae =2{PgP, + e‘^^QgQr), 

A. = [P2 + e-2-P| + e^P{Ql + e-^^Qj)]. 

Notice that the wave equations evolve in an unconstraint manner with respect to (^. 

These equations serve two purposes. They are used to generate initial data, and can 
be readily integrated for obtaining reference solutions against which the simulations can be 
compared. A simplectic integration scheme has been given by [Q . For reasons of convenience, 
a leap frog integration scheme in combination with spectrally accurate spatial differentiation 
by the Fast Fourier Transform has been used here. Thus, general Gowdy waves are at hand 
with high accuracy as a reference in studying the simulations in the tetrad approach. In the 
polarized case, furthermore, analytical solutions are available. 


A. A polarized Gowdy wave 

Polarized Gowdy wave result in the special case Q = 0, so that the Gowdy line-element 
becomes 

= e~ (—e“^’'dr^ dd^) e“^[e'^dcT^ e“^dd^]. (6) 

The reduced wave equation is linear in P, and can be solved by the method of separation 
of variables, giving 

P(r, 9) = (a„ cos n9 + sin n9 ), (7) 

where Zo is a solution to Bessel’s equation of order zero, and where the linearly growing 
solution P = Pq -|- rPi (Pq and Pi constant) has been suppressed. 
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A polarized Gowdy which evolves towards a singularity is exemplihed by (Equation (4.2) 


in 


i) 


Po{T,e) = Yo{e cose, 


( 8 ) 


where Iq is the Bessel function of the second kind of order zero. With <5 = 0 the Hamiltonian 
constraints for A at r = 0 become [Q 


so that (H) implies 


— ‘^PbPt, 

Xr = P^ + e-^^Pl 


(9) 


X{t, 0) = ^Yo{e '^)Yi{e ^)e ^cos26* + r(r), 

T{r) =y^-.{Y^{t) + Y,^{t)}tdt. 

Equations (||) and provide for initial data for simulations, an are an analytical reference 
solution for polarized waves. The computational procedure to treat the integral T(t) is 
outlined in the Appendix. 


B. Unpolarized Gowdy waves 


Unpolarized Gowdy waves follow from general initial conditions on the P{t, 9) and 
Q{t,9). Following the earlier computations of Berger et el. p], we consider 


p(0, 6) = 0, drP{T, 9) = A cos{9), 

Q{d,9) =Bcos{9), drQ{T,9) =0, 
where A and B are constants. We then have the expansions 


( 11 ) 


P{t, 9) r.^ATcos9-XT^B^sm^9, 

Q{t,9) ~ (1 — iHr^) COS0, (12) 

A(r, 9) ~ r(A2 cos^ 9 + B^ sin^ 9) - t‘^B^ sin^ 9 

valid to second order in r. These expansions are used to compute initial data for the 
simulations. 
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III. INTEGRATION SCHEME 


The integration scheme combines an evolution scheme for the tetrad connections Ua^u and 
an evolution scheme for the tetrad legs (e^)^. The scheme is one-dimensional, second-order 
accurate in the time-coordinate r and spectrally accurate in the space-coordinate 6. The 
scheme uses a straightforward leap frog time-stepping algorithm, and, in view of periodicity 
in 9, differentiation by the Fast Fourier Transform. In this fashion, the errors are essentially 
due to the time-stepping algorithm. 

In what follows, we shall work with both tensors and their densities, denoted by a tilde, 
e.g., (f) = Time is discretized as = nAr, indicated by a superscript: 0(r„) = 0"'. 

Indices a, b will refer to all four space-time indices, and p, g, r, s, u, v to the spatial 6*, 5 and 
a only. 


A. Evolution of the connections ca, 


afiiy 


The wave equations for Uafiu are implemented directly through the divergence of the Rie- 
mann tensor, its representation in terms of the connections and the Lorentz gauge condition 


Cfiu V 0- 


The integration is based on the four-divergence of the Riemann tensor in the tetrad ap¬ 
proach, = 0. In coordinate notation, this gives = 0. Application 

of leap frog time-stepping gives the iterations 


= {R^'' 


fll// 


(13) 


-2AT{d,W^, + K, 

Evolution equations for the connections Up^i, (where p is a spatial coordinate) then follow 
from the identity RrpfMu = dT^piJiv — dp^TtJ,u+[^T, ^p]piv Application of leap frog time-stepping 
gives 


[u. 


Pliu) 


^n+1 _ 


= o;. 


\n—1 


Pliu) 


“t“2 At(Rt-p^^ “1“ dpUJ [UJr 5 ^p\ Ifiiy') 


(14) 
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Evolution of the remaining are determined from the Lorentz gauge conditions = 
— O’ hence 

(15) 

-2ATa,(i9>;„)“ 

Note that appears in Step (b) of the iteration scheme. An update is alge¬ 
braically related to and {Rpq^u)"'~^^, the hrst of which follows from Step (a) and 

the second of which follows from Step (b). Suppressing momentarily the tetrad idices fiu 
and the time-label n -|- 1, we have 

Rrp QrcddpR ‘^9T[T9r]pR R 9T[s9r]pR 

= 2gr[r9r]pR'"'' + 9T[s9r]p9^°‘9^''Rab ( 16 ) 

= ApR-rv + Bp, 

were 


Al = 2gr[,gr]pg^^g^-, 

Bp = 2g^i^gr]pR'"^ + 

Note that Ap and the second term in Bp is due to a shift, grp] in particular, Ap 
grp = 0. It follows that 


( 17 ) 

0 whenever 


{6; - A;)Rr, = Bp. (18) 

As a system of 3x3 equations, this is readily inverted. With all quantities at tn+i = (n-|-l) At, 
(^) determines updates {Rtpij,u)^~^^ from the results of the previous steps. 


B. Evolution of the tetrad legs (e^)^ 

The tetrad legs satisfy the equations of structure 

d[a{eg)b] = {^u)[b^a]p ■ (19) 

Clearly, these equations leave the tetrad lapse functions = (e^),- as free variables. They 
are related to the Hamiltonian lapse and shift functions {N, Np) through 
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N^(e‘‘),= (N^N’’-N^,N„). 


( 20 ) 


dar 

These tetrad lapse functions govern the evolution of the tetrad legs in the equations of 
structure: 


dr{ef,)b + uJrll = dbNf, + = dbN^. 


( 21 ) 


Thus, conditions on Qar from the line-element (|^) result in a system of implicit equations 
for the lapse functions N^, as discussed in Section IV. 

Integration of (^) can be written using the fundamental matrix A^{t]To): a hnite 
Lorentz transformation satisfying 


I A;^(r,r) =h;(, 


with which the solution to (^) is 


( 22 ) 


{ef,)b{T) = A^{T-,To){e^)b{To) 

+ IroK(^'^ s)dbN^{s)ds. 

Here, only the r—dependence has been made explicit. Indeed, if the convolution integral on 
the right hand-side were to vanish, the evolution of the tetrad becomes a pure SO(3,l,R)- 
gauge transformation, leaving the metric gabir) = (eQ,)a(r)(e")b(r) constant. The general 
representation (^31) shows explicitly that the evolution of the metric takes place in the 
(flat) tangent bundle of the four-dimensional, physical manifold. The integral in (p3D is 
implemented using Gaussian integration, including a Taylor expansion for the integrand as 
provided by the updates in Section 3.1. Thus, in the case of given connections, the tetrad 
evolution is carried out with machine precision. In general, a Taylor expansion of up 
to its n—th derivatives at tn in ( P3|) provides an update of the tetrad legs with n + 1—th 
order accuracy. 
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C. Equations for A^{t]Tq) 


In the Gowdy line-element, Fg^ and R^-oc have a 2x2 block diagonal structure 


as matrices in c, d. Initial data for the tetrad (Section D) can be chosen with (e^)b likewise, 
when regarded as a matrix in yu, 6. This obtains tuefiu, and Rrefiu in block 

diagonal form at r = 0. Consider now the evolution problem for the connections in terms 
of drU} 0 f,^ = deUr^u + RTe^lv - and the Lorentz gauge, = 0. Here, R^e^u = 

RTecd{.^ 0 .Y{ei^y, where (the structure of) Rrecd is as mentioned above, and where the tetrad 
elements satisfy ( PD with = (*,*,0,0) (Section IV). Consistent evolution is obtained in 
which the block diagonal structure of the forementioned tensor elements is preserved with 


aa 0 


(24) 



where a = a{T, 6 ) and b = b{T, 9), and 



^0 ( 0 - 1 ^ 


By explicit exponentiation, the Lorentz transformation of boosts on the {ctY and (ee)^, and 


rotations on (es)^ and (ca)^ is therefore 



= 


(26) 


Here, we adapt the nomenclature of 0, with the 2x2 matrices K and L given by 



(27) 


with 



(28) 
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D. Initial data and choice of tetrad 


The integration scheme requires uja^iv and drUJafj.u for obtaining initial values of the Rie- 
mann tensor on r = 0. For the polarized Gowdy wave (Q = 0), we can simply use 

(eT)‘ = (ei0,0,0), 

(ee)‘ =(0,ei,0,0), 

(29 ) 

(es)' =(0,0,e-^,0), 

(eA)' =(0,0,0,e^). 

In case of the full Gowdy metric {P,Q ^ 0), care must be given to the choice of ini¬ 
tial tetrad to ensure sufficient differentiability for obtaining smooth connections and 
Riemann tensor Rabfiu- The tetrad legs for the t — 9 coordinate plane can be chosen as 
before, 


(ct)^ — (e4,0, 0, 0), 
(ee)' = (0,eT0,0). 

The tetrad legs for the a — 6 coordinate plane, 

(cs)cr { 6 a )(7 

h/ = 

(eA)iy 


( 30 ) 

( 31 ) 


( 32 ) 


must satisfy EE~ — G, EG ^E~ — 7, where 7 is the 2x2 identity matrix and G is the 
matrix with the metric components 


1 Q ^ 

yq Q^ + e-2^ ^ 

E can be defined uniquely as the symmetric positive definite square root (c/. 
this end, let 



( \ 

/ 

G = 

Qaa 9 <jS 

= 


^ QaS 9ss y 

V 


(33) 


ofG. To 


m 


/ 

V 


cos (p — sin (p 
sin p cos p ^ 


( 34 ) 
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denote both the matrix containing the tangents to the tetrad legs, and the rotation matrix 
on the tetrad indices fi = a,6. Then 


E = A{4>) 


( i \ 

' \l O' 

0 A! 




(35) 


where the eigenvalues A± of G can be given in terms oi z = ^{Q^+e ^^—1) and r = 
as 


A-t = (1 + 2 ; ± r)e 


-T + P 


( 36 ) 


The rotation angle (f) follows from the equations 


2 ; sin 2(j) + Q cos 20 = 0, z + r cos 20 = 0. 


( 37 ) 


The second equation in (|37|) gives cos 20 = — and hence sin20 = Without loss of 
generality, cos0 > 0, so that O<0<|if(5>O, and —|<0<Oif(5<O. The additional 
rotation y4(—0) on the tetrad indices in the right hand-side of (^) effectively regularizes 
E as G approaches the 2x2 identity matrix (when P and Q become small, in which case 0 
becomes ill-dehned). 

The connections oja^xu are now computed using its metric dehnition 


(c^)c^a(6j/) -|- 


(38) 


Initial data for the Riemann tensor of Section Ilia now follow. 


IV. GAUGE CONDITIONS FOR GOWDY WAVE 

The Gowdy line-element prescribes a certain slicing of space-time, which are incorporated 
in the simulations to enable comparison with the analytical solution for the polarized and 
pseudo-unpolarized Gowdy waves, and the reference solutions obtained by integration of the 
reduced wave equations (^ and Hamiltonian constraints (^. 

The Gowdy line-element (j^) has algebraic slicing conditions 
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9 tt C 909 ^ 

9t9 0) 

9 Ta 0 , 

9t8 0 - 


We have 


( 39 ) 


9rT = (ea)r(e")r = - + A^e, 

999 = (ea)0(e")6) = —(er)i + (69)0, 

9 t 9 = (ea)r(e")6) = — A^T (er)6» + (69)5). 

The shift condition gr 9 = 0 gives 

Nq _ {eT)e 
Nt (ee)^’ 

so that 

m = (ee)? (l - |g]|) = (ee)i (l “ ■ 

Becanse also 


( 40 ) 


( 41 ) 


( 42 ) 


9 tt 



the lapse condition grr = —e g 99 becomes 




{ee)l 


(43) 


(44) 


We thns have the following conditions on the tetrad lapse fnnctions 

/ 

Nt = -e-'^{ee)9 = -e-^a^{e^)0, 

Ne = -e-^{eT)9 = -e-^a^{e^)e, 

(45) 

= 0 , 

A^a = 0 . 

Upon snbstitntion of (^5[ ) into the eqnations of strnctnre (P5[), two implicit eqnations for the 
( 69)0 and (e'r)^ resnlt. 
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A. Solution of gauge conditions 


The nontrivial t — 6 gauge conditions in (^31) are separated from the trivial gauge con¬ 
ditions on the a — 6 coordinates. Our starting point, therefore, is 


{e„)eiT,0) = K"{T-,Ta,0)r]^{T,e), {p,v^ T,e), 


(46) 


and to work in the t — 6 sector only. In what follows, Greek tetrad indices run through T, 0. 
Also, the spatial (0—)dependence will made explicit only when needed. 

The original gauge conditions (^) are given in terms of the metric, which are invariants 
under pointwise K transformations applied simultaneously to both (ct)^ and (ee)^. Upon 
substitution of in (45), therefore, K can be factored out, leaving a linear equation for 
Vt.- 

Expressed in (^61), the gauge conditions (|45|) become 




(47) 


Substitution of (^61) and (^) into the equations of structure (^) gives two implicit equations 
for r]fj, = {r]T,rie). Note that 


(48) 


=(3{T,e)a^ 

dBK^{T- To, 6) = a^K^^{T- tq, 0)a(r; tq, 9), 

where (5 = P{t,9) is a coefficient function, and Q;(r;ro,p) = d 0 a{s,p)ds. Together with 
dtK^ = —aa^K^ from the second equation in (p^, and the algebraic properties = 

A^ajj' and cr^a^ = h)), the linear equation for 77 ^ follows: 


drVf, + e c(r, To)r]^) = 0, 


(49) 


where c(r, Tq) := /3{t) -f Q:(r; Tq). Because a^r; r) = 0, c(r, Tq) combines with dgu^^ to 
hrst and second order in (r — tq), respectively. 

Momentarily suppressing the tetrad index p, we now look for a solution to (|49|) in 


^2. 1 


V = Vo + {r- To)rii + ^(^ - 'roTv2 + “ ^o)^h3 + 


(50) 
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where each rjk = 'rik{9) has only 0—dependence. Similarly, we write 


gr ^ gro + (7- _ ro)e^° + |i(r - ro)^e^° H- , 

c(r, To) = Co + (r - ro)ci(ro) + (r - ro)^C 2 (ro) H-. 

Expansions (|^) and (^T]) can substituted into (^9]); matching coefficients gives 


( 51 ) 


+ Coho = 0 


(52) 


for A: = 0, and 

/c 2 k—1 

H uVk+i-i + T(^^deri^k + Cfcho + t — icm-i = 0, 

1=0 ^ z=o ^ ^ 

for k > 1, i.e., 

Vk+i+ Ef=i jiVk+i-i + e~^°{l(^ii^der]^k 

+CkVo + ]^iCiVk-i) = 0. 

The hrst few terms are 

/ 

Vk-o = (e/^)e('ro), 

h^i = -e~'"°{a^de'r]uo + coh/.o), 

hM2 = -h^i - e~'^°{a^deriui + cit/^o + coh^i): 

h/iS Vk'^ 

-e"^°(ia^'"(9e77j.2 + C2hM0 + cih^i + |coh/.2), 


(53) 


(54) 


(55) 


This series expansion has been implemented numerically using a numerical cut-off of l.D-12, 
thereby maintaining the Gowdy line-element within machine accuracy. 


V. RESULTS AND CONCLUSIONS 

The performance of the numerical implementation has been studied by varying both 
space and time discretization. Verihcation of second-order accuracy has been obtained by 
computing the ratio of the errors (as a function of time) of those with 2 m time-steps to those 
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with m time-steps for a given final value of r. In all computations, a moderate degree of 
discretization in the 0—variable has been found adequate, because of the spectral accuracy 
in 0—differentiation by FFT. The computations have been performed in the Gowdy line- 
element, using the freedom of space-time slicing in the present formulation. Numerical 
results and convergence data have been obtained for the polarized, pseudo-unpolarized and 
the two unpolarized Gowdy waves. The results are shown in Figs. 1-3, for both the polarized 
wave (Figure 1) and two unpolarized waves (Figure 2,3). The results show second-order 
convergence in accord with the numerical scheme. Also shown is convergence in the errors in 
the Ricci tensor (also second-order). The results suggest that more advanced time-stepping 
algorithms should be applicable, to adapt for for more accurate, long-time computations. 
At this stage, the results do not indicate a need for implicit time-stepping. 

The computational problem of binary coalescence of neutron stars or black holes requires 
a continuing development, including an extension to three-dimensions, adaptation to non¬ 
periodic boundary conditions, and possibly a minimization of the error in the Ricci tensor. 

Two directions for further exploration stand out, which are motivated by the structure of 
the equations. Firstly, the present formulation offers the hrst strictly-hyperbolic formulation 
with complete freedom of foliation (in particular, there is no restriction on the equivalent 
Hamiltonian lapse function). It will be of interest to exploit this for new avenues in the 
treatment of horizon boundary conditions, and possibly in the extraction of wave forms at 
the outer boundary as well. Secondly, the implementation is based on the four-divergence 
of the Riemann tensor and contains no Ghristoffel symbols. Somewhat analogous equations 
are given by Faraday’s equations, which recently have been successfully implemented using 
cylindrical coordinates in simulations of magnetized relativistic jets P,|^. Here, the noto¬ 
rious axis instabilities are fully regularized P, whence it will be of interest to extend this 
regularization to the present equations. Gylindrical coordinates would be desirable in view 
of the expected simplihcation in the asymptotic wave form at large distances. 

Finally, inclusion of a non-trivial stress-energy tensor is required in the problem of coa¬ 
lescence of neutron stars and possibly black holes as well 
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Appendix. The analytical solution of the polarized Gowdy wave is given by Po{t,9) in 
(I) and A(r, 6 ) in (0). It has been found useful in debugging the program and in studying 
convergence to be able to monitor each and every variable during the simulation in compar¬ 
ison with the exact solution. To obtain these and all derived quantities, such as uJa^iu and 
Ratin', at every time-step, numerical evaluation of the right hand-sides of (§) and (|^ are 
needed. A convenient method to do so is by a direct evaluation of their dehning expressions 
in terms of the metric, using numerical differentiation by hnite differences across small time- 
and space-intervals e (much smaller than the time-step size and spatial discretization in the 
numerical integration of the wave equations). This requires accurate evaluation of 


nr) =yn{y^n+yo\t)}tdt 

= Uo{n\e-n+yne-n}e-^^dr 


(56) 


up to its second derivatives about each of the = nAr. 

Generally, T(r) on [tq, Tb] is given accurately up to its k-th derivatives in each subinterval 
Dn = [y — ke, Tn + ke] about = nAr = by a chain of function elements {fn, Dn} is 
used (see, e.g. [TO). Here, the fn are polynomial approximations 


/„(t) ^ CnO + C„i(t - T„) + Cn'iir “ H- 


(57) 


to T(r) on Dn, and the Dn cover [—ke, mAr -|- ke]. Naturally, the degree of the fn is chosen 
to be greater than or equal to k. In the case at hand, k = 2, and the degree of fn is chosen 
to be three. 
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The coefficients Cnu ^ = 1, 2, • • ■ follow by direct evaluation of the (analytic) prescription 
of the integrand of T{t) and its (analytical) derivatives. A set of accurate values of the c„o 
then follows by hrst evaluating {T(r^)} on 


/ ^ A / / .1 \ Stt 

= -mAr (1 + cosxfc), Xk = —, 


(58) 


using integration of its integrand on the {r(,} by the Fast Fourier Transform, where N is 
the degree of the FFT, followed by an interpolation of the T(r^) to the T(rn) using the Cni, 

In our computations N = 256, which gives results on T(r) to within machine accuracy. 
Of course, for every derived quantity involving numerical differentiation, there is a loss of 
accuracy determined by the numerical differentiation parameter e. We have chosen e = 10“®. 
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FIGURES 

FIG. 1. Shown is the simulation on 0 < r < 5.12 of the polarized Gowdy wave. Distributions 
X{t,6) and P{t,9) are displayed (upper windows), together with their 0—distribution at r = 5.12 
(middle windows). Errors are obtained in computations with a consecutive doubling of the number 
of time-steps {m = 512,1024, 2048), and are given as a function of time in the lower window to the 
left for both P and A (o for m = 512, x for m = 1024 and lines for m = 2048). The errors show 
proper second order convergence (lower window to the right), obtained upon dividing the errors 
for m = 1024 by those of m = 512 (x), and those for m = 2048 by those of m = 1024 (lines). 
These errors have been computed with reference to the analytical solution to Gowdy’s reduced 
wave equations. 

FIG. 2. Shown is the simulation on 0 < r < 5.12 of the unpolarized Gowdy wave in response to 
initial data with A = B = 1. Discretization in 0 is n = 64 (32 points displayed). Errors are obtained 
in computations with a consecutive doubling of the number of time-steps (m = 512,1024,2048), 
and are given as a function of time in the lower window to the left for each of the three functions 
P, Q and A (o for m = 512, x for m = 1024 and lines for m = 2048). The errors show proper 
second order convergence (lower window in the middle), obtained upon dividing the errors for 
m = 1024 by those of m = 512 (x), and those for m = 2048 by those of m = 1024 (lines). The 
errors in the Ricci tensor for the three discretizations (lower window to the right) likewise show 
second order convergence. An additional dashed line further shows the error in the Ricci tensor for 
extremely high time-discretization, m = 32768. These errors are computed using numerical results 
from Gowdy’s reduced wave equations for comparison. 

FIG. 3. Shown is the simulation on 0 < r < 7.68 of the unpolarized Gowdy wave in response 
to initial data with A = 0 and B = 1. Discretization in 0 is n = 64 (32 points displayed). The 
windows display the same variables and errors as described in the previous figure, with errors again 
computed using results from Gowdy’s reduced wave equations for comparison. A notable transition 
occurs in the behavior of the errors, specifically in those of the Ricci tensor, as steep gradients are 
formed in P and Q (and to a lesser extend in A as well) about r = 5. 
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